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We present finite temperature analysis of a quasi2D dipolar gas. To do this, we use the Hartree 
Fock Bogoliubov method within the Popov approximation. This formalism is a set of non-local 
equations containing the dipole-dipole interaction and the condensate and thermal correlation func- 
tions, which are solved self-consistently. We detail the numerical method used to implement the 
scheme. We present density profiles for a finite temperature dipolar gas in quasi2D, and compare 
these results to a gas with zero-range interactions. Additionally, we analyze the excitation spectrum 
and study the impact of the thermal exchange. 
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I. INTRODUCTION 

Two hot topics in ultracold physics are two dimen- 
sional (2D) gases and dipolar gases. Reduced dimen- 
sionality enhances the quantum mechanical character of 
the system. In a homogeneous 2D system, quantum 
phase fluctuations are so strong that at finite tempera- 
ture phase coherence cannot be established and conden- 
sation does not occur. There is, however, an intrigu- 
ing phase transition called the Berezinskii-Kostcrlitz- 
Thouless (BKT) transition [1] which occurs when the 
temperature is lowered and there is no longer enough 
thermal energy to unbind vortex and anti-vortex pairs. 
This binding of vortices reduces phase fluctuations so 
that quasi-long range order can form [2]. Interestingly, 
a trapped 2D gas can Bose Einstein condense (BEC) at 
finite temperature [2, 3]. The BKT and BEC transi- 
tion have been studied in trapped ultracold gases, first 
by observations of phase defects [4] , direct observation of 
vortices [5], and then changes in the density profile due to 
the onset of superfluidity [6] . Additionally, universality 
has been observed near the BKT transition [7] . Success- 
ful methods which study such systems are Monte Carlo 
[8, 9], mean field [11, 12] and classical-field methods [13- 
15]. 

In dipolar systems the interactions are non-local and 
the possibility of creating strongly correlated gases are 
tantalizing, especially in quasi-2D (q2D) where zero point 
motion in the axial direction is included. For example, 
some studies have used Monte Carlo methods to study 
phase transitions such as crystallization [16]. This study 
included temperature as a variable and observed melting 
of the dipolar crystal. The BKT transition of a homo- 
geneous dipolar gas has been studied with Monte Carlo 
methods to examine superfluid fraction and excitation 
spectra [17]. Other studies have looked at the structure 
of the dipolar gas and the impact of temperature on the 
phase of the gas [18, 19]. Such studies have focused on 
strongly interacting dipolar systems. However, the bulk 
of the work on q2D gases has been at zero temperature. 
For example, phonon instability [20] and anisotropic su- 
perfluidity [21] have been predicted. 

Recently, chromium (Cr) and dysprosium (Dy) have 



been Bose condensed and exhibited strong dipolar effects 
[22, 23]. Additionally, progress towards a q2D dipolar 
gas has been made with a layered dipolar system for Cr 
atoms in a one dimensional optical lattice [24]. With 
such progress on dipolar gases, direct study of q2D dipo- 
lar systems has begun. An important question is how 
will dipolar interactions impact the quantum behavior 
of a q2D gas? The long range nature of the dipolar in- 
teraction may lead to interesting physics for the thermal 
system, especially relating the BKT transition and phase 
coherence. Additionally, little work studying the impact 
of temperature on q2D dipolar gases has been conducted 
for reasonable experimental parameters. 

With an eye toward this unexplored physics and aiding 
experiments, we study the finite temperature physics of 
q2D trapped dipolar gases. We use the standard beyond 
mean-field method: the Hartree Fock Bogoliubov within 
the Popov approximation (HFBP) [25]. The HFBP has 
been successful in studying ultracold atoms, for example 
it has been used to explain the temperature dependence 
of collective excitations in a 3D BEC [26]. In the case 
of 3D dipolar gases, Rcf. [27] used the HFBP to study 
temperature effects on the biconcave structure dipolar 
gases can have [28]. This study neglected the thermal 
exchange. Other recent work used mean field methods 
to study the stability of finite temperature dipolar gases 
[29]. 

The aim of this paper is two fold. First, we develop 
a numerical method to solve the non-local HFBP. This 
method includes the non-local interaction and exchange 
effects. We show in detail how the interaction and ther- 
mal exchange effects arc included; it relies on a parallel 
implementation of the method. The method can also 
be applied to dipolar fermions. Second, we compare the 
contact gas and the dipolar gas with the HFBP. In the 
comparison of the two gases, we look at the condensate 
number as a function of temperature. We also articulate 
the role of thermal exchange in the gas by solving the sys- 
tem with and without thermal exchange. We compare 
the density profiles of the gas at various temperatures. 
We look at the excitation spectrum as a function of both 
temperature and total particle number. We also classify 
the lowest excitation modes. To the author's knowledge, 
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these results are the first for a q2D trapped-finite tem- 
perature dipolar gas with thermal exchange effects. 

II. EQUATIONS OF MOTION 

We study a quasi-2D dipolar system at finite tempera- 
ture. To do this, we will employ the Hartree Fock Bogoli- 
ubov method within the Popov approximation with non- 
local interactions [25, 30]. The HFB breaks the second 
quantized bosonic field operator into a condensate and 
non-condensate (thermal) component: ^ = [^/Wo(f) {x) + 



9{x)] where we have replaced So — > \/N a . Here x repre- 
sents all required coordinates, for the 2D case it will be 
p. We will use the Bogoliubov transformation for the 
thermal part: 6(x) ~ J2 a i u »( x )^ae^' lu ' a ' t ~ v^(x)d* a e lulat ] 
where d (d*) is the bosonic annihilation (creation) oper- 
ator for a quasi-particle (hole). They obey the bosonic 
commutation relations: [a a ,aV\ = <5 Q( g and [a,,,^] = 
[d*,dg] = 0. To derive the equations of motion we start 
with the Heisenberg equation of motion with a non-local 
interaction [30]. We find for our system the modified 
non-local Gross Pitaevskii equation (GPE) is 



H<p (x) = (h + J dx'n(x')V^j <f) (x) + J dx'lhix,^) - m(x,x')]V(j> (x'), (1) 

where Hq is the kinetic energy and trapping potential. The interaction potential, V(x — x'), has been written as V 
for simplicity. The total density, n(x) = uq(x) + fi{x) is made up of condensate density and thermal density. The 
non-local correlation functions are: no(x,x') = N <f)Q(x)(f)o{x') and n (x) = n (x 7 x); n(x,x') = (9* (x)9(x')) is the 
thermal correlation function and the local thermal density is h{x) = hix^x). rh{x,x') = (9(x)9(x')} is the anomalous 
thermal correlation function, and in the Popov approximation is neglected. The quasi-particle wavefunctions are 
required to construct the thermal and anomalous correlation functions. At equilibrium, these wavefunctions can be 
found from a non-local Bogoliubov-de Gennes equation: 



+ J dx'm(x, x')Vv a {x') — J dx'm (x,x')Vv a (x') 



(2) 



-u n v n 



+ 



J dx'rh*(x, x')Vu a (x') — J dx'm^(x,x')Vu a (x'). 



The total number of atoms is N = N + N and N — 
J dxn(x). The normalization of the quasi-particle wave- 
function is 1 = J dx(\u a \ 2 — \v a \ 2 ). The anomalous con- 



densate correlation function is m (x,x') = N <j)(x)(j)(x'). 
The thermal correlation functions are defined in terms of 
the quasi-particles: 



h(x,x') = J2K( X > a ( X ) + v a (x')v* a (x)}N^ + v a (x')v* a (x) (3) 

a 

m(x,x') =J2K(x')v* a (x)+v* a (x')u a (x)]N? e + u a (x')v* a (x), (4) 

i 

where Ng e = (Ze hUa / kT — 1) _1 = (a* d a ) where k is cation is made. 
Boltzmann's constant, T is the temperature, and Z = 
1 + 1 /No . This relation is simply the Bosonic occupation 
of thermal modes. These equations are consistent with 
Refs. [27], [30], and [31] when the appropriate simplifi- 
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To make the discussion of Eqs. (1) and (2) simpler, we 
introduce a compact notation: 



is: 



{H + D + X + M)fa= vufo, 
(H Q - p, + D + X)u a + Mv a 
(H - p. + D + X)v a + Mu a 



(5) 
(6) 



-i0 a V o 



The first term is the standard kinetic and potential 
energy, Hq. For the remaining interaction terms, we 
use the capital letters to indicate the inclusion of the 
x' integral, so the direct interaction is D = Dq + D 
(condensate and thermal), and for example D (x) = 
J dx'[no(x')]V(x — x'). The exchange interaction is 
X = Xq + X (condensate and thermal), and for example 
X(x,x') — J dx'[h(x,x')]V(x — x'). Finally, the anoma- 
lous correlation function is M = M—M , and for example 
M(x,x')= J dx'[rh(x,x') - m (x,x')]V(x - x'). We will 
use the Popov approximation, which assumes M = 0, 
leading to M — — A if (f>o is real. This approximation 
keeps the spectrum gapless. 

We project out the condensate mode from the quasi- 
particles, as was done in Refs. [32, 33] which uses 
Q = 1 — |0)(0| and (a;|0) = (j>o(x). The projection op- 
erator is applied to all operators that are not in H GP = 
Hq — /i + D + X, such as Mq and Xq, and for exam- 
ple, Mq — > QM Q. This projection method keeps the 
spectrum gapless. 



III. THE INTERACTION 

We are interested in the q2D dipolar system assuming 
no axis of symmetry, so we take the dipole moment to 
be d = d(zcos(a) + xsm(a)), and we assume harmonic 
confinement and that ui z 3> lo p where the uii are the 
trapping frequencies ( u> x = ui y = u> p ). This allows us 
to assume only one transverse mode is occupied in the 
z direction. We can then evaluate the interaction by 
factoring the wavefunction as ij)(r) — xo(z)ip(p), where 
tp is either (f> , u a , or u a and xo is eigenstate of H in 
the z direction. We are then able to integrate out z, and 
obtain an effective interaction. We will also consider a 
contact interaction. The interaction we use is: 



V(p)=gS(p)+g d V dd (p), 



(7) 



where g = \Z8nh 2 a s /ml z is the strength of the contact 
interaction, a s is s-wave the scattering length (a s -C l z ), 
l z = yjh/muj z is the axial harmonic oscillator length, 
gd = d 2 j\phxl z and d is the dipole moment of the par- 
ticles. This interaction leads to the dipolar length scale: 
Id = md 2 /h 2 . In practice, we rescale the equations 
of motion into oscillator units, where t he unit s of en- 
ergy and length are: Hco p and l p = yjh/muj p . After 
this rescaling, the dipolar interaction strength becomes 
ga = ld/{lzV%n)- Finally, the dipolar interaction in q2D 



V dd (p-p)= J dzxl(z) J dz' X 2 (z')V! D (r-n%) 

V2 d D (r)=d 2 (l-3(d.ff)/r*. 

The interaction is evaluated in momentum space. As an 
example, the direct interaction is: 



D(p) = J- 1 \n(k)V (fc) 

T~r/T\ 47T „ / kl z \ 

Vik) = 9d Y F [f2)> 



(9) 



where T is the Fourier transform operator and n(k) = 
F[n(p)\. V k is the momentum representation of the in- 
teraction. The function F(k) is the fc-space dipolar in- 
teraction for the q2D geometry It has two contribu- 
tions coming from polarization perpendicular to and in 
the direction of the dipole tilt, F(k) — cos 2 (a)F±(k) + 
sin 2 (a)F||(fc) where a is the angle between z and the 
polarization vector d. Its contributions arc F\\(k) = 

— 1 + 3v / 7r(fc 2 /fc)e fe erfc(fc), where k d is the wave vec- 
tor along the polarization direction in the x-y plane, 
erfc is the complementary error function and F±(k) = 
2 - 3V^ce fe2 erfc(fc) [20, 21, 34]. 

To evaluate the exchange interaction, we use a simple, 
yet memory intensive, solution: we explicitly construct 
the interaction on the same space-space grid as the cor- 
relation functions. We put n(p,(7) on a grid specified by 
the indices s and t, representing p and p . Each grid has 
n s spatial points, and a correlation function has n 2 grid 
points. 

To evaluate the interaction, we work in momentum 
space. The interaction on the space-space grid is: 

V st = Wj k V k W kt . (10) 

where W ks is the operator that transforms between space 
(s) and k-space (k) representation via the spectral ba- 
sis set [35] . The operator takes the place of the Fourier 
operator in our algorithm. An important aspect of the 
method is that it regularizes the interaction and avoids 
the logarithmic divergence as p — > ff encountered by di- 
rectly evaluating Eq. (8) because W kt involves a project- 
ing onto the Gauss-Hermite basis set [35]. Constructing 
the interaction involves the most costly step, requiring 
rig operations per node. This is performed only once for 
a set of interaction parameters (l z and a) and basis set 
size. 

We have constructed a parallel implementation of the 
method. We distribute the s index (p grid) across the 
computing nodes (p); we will denote this distributed in- 
dex as s p . This way, each core handles only a fraction 
of the correlation functions and interaction. In the fu- 
ture, we will construct the interaction in the partial wave 
expansion: V(p) = ^2 m V m (p 7 (f>). This expansion will 
be necessary to handle an interaction which depends on 
quantum numbers such as n z (confinement in z) or spin. 
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IV. NUMERICAL IMPLEMENTATION 

We solve Eqs. (5) and (6) by expanding the wavefunc- 
tions on the basis that diagonalizes Hq. We will focus 
on the harmonically trapped case where H^Xa = ^aXai 
I dpx* a (p)Xp(p) = <W, and e a = fuv p (m x + m y + 1) 
where m a is an integer. In a basis set, the evaluation 
of the direct interaction term, D a p, is straight forward: 
I dpXa(p)[J dx'V(p- p , )n(p , )]xp(p), where the quantity 
in the square brackets is the effective potential from the 
interaction. The more complicated terms are the non- 
local exchange terms, such as the thermal exchange: X a p 
which is / dpx* a (p) J df?n(p, ff )V(p- ?) X p{P). 

The non-local exchange term in the spectral basis is: 
X aa = U^ s V st n st Uta, where U SIJ is the transformation 
between the spectral basis (a) and space basis (s). The 
numerical procedure to evaluate this on node p is: (1) 
Y Sp t = V Sp tn Spt , multiply the interaction and correlation 
function; (2) h Sp<7 — Y Spt U tcr , project the t spatial ba- 
sis onto the a spectral basis; (3) X^ a = J2 Sp UT p a^s p a, 
project s p spatial basis onto the a spectral basis; and (4) 
Xao- = Sp^ao-i collect and sum. In comparison, the 
direct terms are straight forward and evaluated in mo- 
mentum space at each step. $ s = W s kVkWktnt, where 
& s is the effective potential. Then we must project onto 
the basis: D aa = £/J s 4> s £7 s(7 , and this is done in parallel. 

The procedure to find the full solution is: set temper- 
ature, T, and iVo and pick a targeted value for the total 
number of particles in the system, N tar get, and set h = 
and m — 0: 

1. Solve Eq. (1) for <p and p (in basis set). 

2. Construct condensate exchange term, X , with 
n n (p, /f) and put in basis set. 

3. Solve Eq. (2) for u a and v a . 

4. Construct n(p, ff) and m(p,p'). 

5. Use the semi-classical LDA to supplement the ther- 
mal tail (see below), which gives N and therefore 
N = N + N. 

6. Construct thermal and anomalous exchange terms, 
X and M (if not using Popov approximation), and 
put in basis set representation. 

7. Adjust A to get the desired N tar get of atoms. We 
require N to be within 1% of A 'target ■ 

8. Go back to step 1 until self-consistency is reached. 
We converge the number of thermal atoms, so that 
\hi — fij_i| < 5 • 10~ 5 , when we are near enough 

^target • 

With the numerical method in hand, we are ready to 
proceed to the physical examples and the results. 

The semi-classical, local density approximation (LDA) 
is very important for large p and high momentum states 
on the grid. At high temperature with a manageable 



grid size (much bigger than the condensate), we find that 
there can be a appreciable number of atoms outside of 
quadrature grid. For the large p contributions, the den- 
sity is low enough that the analytic, non-interacting so- 
lution can be used to account for these off-grid particles. 
We integrate the analytic solution over a large region out- 
side of the grid. We have tested that this does not impact 
the convergence of the end result. We include the impact 
of higher transverse modes in the same analytic manner 
on and off the grid. We also assume that these particles 
in the higher z modes are non-interacting. This is rea- 
sonable, especially because symmetry prevents the first 
excited transverse mode from colliding with the conden- 
sate with ui — symmetry. Additionally, a noticable 
number of thermal particles are in momentum states be- 
yond calculated quantum spectrum [36, 37]. To include 
these thermal atoms we use the method from Ref. [27] . 



V. PHYSICAL EXAMPLES 

We will study Cr and Dy atoms. These atoms have 
been Bose condensed and shown to exhibit strong dipolar 
effects [22, 23]. The first example we choose is: 2 • 10 3 Cr 
atoms with Id = 45a where ao is the Bohr radius. This 
gives gd = 0.0028 and g = in trap units. Addition- 
ally, we pick l z /l p = 0.1, this is from using (oj p ,uj z )=2ir 
(16,1600) Hz. We have found that our numerical proce- 
dure using 465 basis states with an energy cut of 30hto p 
converges. To compare this to a contact gas, we match 
the chemical potential at zero temperature by varying g. 
We will only consider a = here, that is perpendicular 
polarization where the dipolar interaction is isotropic and 
repulsive. We also study Dy with a trapping potential of 
(<jJ p ,u z )=2-k (10,1000) Hz. With l d = 400a , this gives 
gd = 0.034 in trap units. This is strongly interacting 
and to maintain the Kohn mode (center of mass slosh 
mode), we have used a low atom number, and for our 
example we use 300 atoms [22]. It is important to note 
that the Cr example is very similar to recent experiments 
[24]. In that work, they were able to vary l z /l p between 
0.1 and 0.17 and have up to 2000 atoms in a layer [24]. 
However, this is a layered system and interlayer dipolar 
interactions are important. 

In contrast to a homogenous system, a trapped ideal 
2D gas can have Bose condensation [3]. In this case, the 
critical temperature is at T c0 /Huj = V&N /n <~ 0.78v^V. 
We will use this temperature to rescale our findings so 
that, to first order, we remove the number dependence 
of the results. This is important because, as we vary 
the temperature, we do not have exactly the same num- 
ber of particles between every calculation. In the ther- 
modynamic limit, the population of the condensate is 
N /N - 1 - (T/T c0 ) 2 [3]. 

For reference, we give several specific numerical exam- 
ples in Table I. For two temperatures, we report p, N /N, 
N, and the Kohn mode energy (should be lfiw p ). For 
the Dy (g d = 0.034) example, at T/T c0 = 0.5, the Kohn 
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TABLE I: Numerical examples of the HFBP for the q2D dipo- 
lar gas and contact gas. 





TIT n 




Nn IN 
iVQ/iV 


N 


rC c\ n n A/f r\c\ p I ni i) 

J.YW11I1 1V1UU.L/ 1 LUJ 


g d =0.0028 


0.05 
0.5 


3.822 
3.457 


0.994 
0.658 


2004 
1994 


0.9992 
0.995 


g d =0.0028 
X — 


0.05 
0.5 


3.814 
3.341 


0.994 
0.660 


2004 
1994 


0.9996 
0.997 


3=0.021 


0.05 
0.5 


3.829 
3.541 


0.993 
0.662 


2004 
1993 


0.9991 
0.9995 


3d =0.034 


0.05 
0.5 


5.124 
4.919 


0.983 
0.525 


300 
298 


0.991 
0.976 


5d=0.034 
X = 


0.05 
0.5 


5.094 
4.524 


0.983 
0.5500 


300 
303 


0.996 
0.987 








FIG. 1: Nq/N as a function of Temperature for the Cr dipolar 
gas (thick black line with circles), the contact gas (green dash 
dot), and for 2000 ideal trapped particles (solid black line). 
Additionally, a gas of 300 Dy particles (dashed red with +) 
and 300 ideal particles are shown (red open squares). 

mode has noticably deviated from 1 by about 2.5%, this 
is about the worst the convergence of this mode gets as a 
function of temperature. For both dipolar examples, we 
show the results for when thermal exchange is negelected 
(X = 0). The lower temperature solutions are similar to 
the GPE solutions at T = 0. 



VI. PROPERTIES OF A Q2D DIPOLAR GAS AT 
FINITE TEMPERATURE 

In this section we study the properties of a q2D dipolar 
gas at finite temperature, and compare it to a contact gas 
and a dipolar system without the thermal exchange. We 
look at several different aspects of these systems. First, 
we look at the condensate number as a function of tem- 
perature. Then, we compare the interaction contribu- 
tions to the total energy. Third, we compare the density 




FIG. 2: The chemical potential for three systems: dipolar 
(black), contact(dashed red), and dipolar without exchange 
(dash dot blue). The interaction terms for each system in 
the Hamiltonian as a function of temperature. The contribut- 
ing terms are the direct condensate ({0o|-Do|</>o)/a*, triangles), 
the direct thermal ({(f)o\D\(f>o) / n, circles), and the thermal ex- 
change {{4> \X /\4>q) / ^l, squares) interaction. 



profiles of the contact and dipolar gas, at various tem- 
peratures. Fourth, we look at the impact of including the 
thermal exchange on the density profile of the gas. Fi- 
nally, we look at the excitation spectrum as a function of 
temperature and total particle number. We also classify 
the lowest excitation modes. 

Figure 1 shows the condensate fraction for the Cr dipo- 
lar gas (black line with circles), the contact gas (green 
dash dot line), for 2000 ideal particles (solid black line), 
300 Dy particles (dashed red), and 300 ideal particles (red 
open squares). The condensate fraction for the dipolar 
gas and contact gas are shifted down from the ideal gas 
(black line), with the dipolar gas being slightly lower than 
the contact gas. For the Dy example, it is clear that the 
interaction strongly depletes the condensate mode. 

In Fig. 2 we look at the chemical potential and its 
contributions from Eq. 1, as a function of temperature. 
The examples shown are dipolar (black), contact (red), 
and dipolar without thermal exchange (blue). We show 
the direct condensate interaction ((0o|-Do|0o)/M' trian- 
gle), the direct thermal interaction (((j>o\D\(j>o) / /i, circle), 
and the thermal exchange interaction ((c/)o\X /\<t>o)/ /i, 
square). For the contact interaction, the thermal ex- 
change is equal to thermal direct interaction, so we only 
show one. For each example, the remaining contribution 
to the chemical potential is from Hq or the potential and 
kinetic energy contributions. 

For the dipolar gas, the thermal exchange term is 
about half as strong the direct contribution. The im- 
portance of the interaction with the thermal particles is 
more important as the temperature is increased. Com- 
paring the dipolar calculation with and without the ther- 
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mal exchange, we see that the chemical potentials agree 
reasonably well, the full calculation has a larger /i at 
high temperature. Looking at the interaction contribu- 
tions, we see that the direct condensate (black triangle) 
and thermal interaction (black circles) make up a smaller 
portion of the chemical potential for the full calculation 
than the one with X = (blue triangle and open circles) . 

In Fig. 3 we show both the total and condensate den- 
sity for (a) a Cr system of 2000 particles and an analogous 
(b) contact system (near equal chemical potentials). We 
then compare the density profiles in (c). The total den- 
sity is in blue (black) for the dipolar (contact) systems, 
and the condensate density is shown in red. The temper- 
atures from top to bottom are T/T c0 — 0.25 (dashed), 
0.55 (dotted), 0.75 (dash-dot), and 0.90 (solid). For both 
the contact and dipolar gas at large p, the thermal den- 
sity behaves as the analytic solution predicts. 

An important point is that the dipolar gas has a higher 
density in the middle of the trap than the contact gas. It 
is hard to see in the figure, but the contact condensate 
atoms have been shifted to the shoulders of the conden- 
sate or near the trap edge. This has implications for the 
temperature at which the superfiuid phase transition will 
occur in a dipolar gas. 

In Fig. 4 we compare both the total and condensate 
density with and without the thermal exchange. We 
look at four temperatures, from top to bottom they are: 
T/T c0 = 0.25 (dashed), 0.55 (dotted), 0.75 (dash-dot), 
and 0.90 (solid). In 4 (a), a Cr gas with (blue) and with- 
out (black) the thermal exchange is shown. Additionally, 
the corresponding condensate density is with (red) and 
without (green) thermal exchange. This example with 
the thermal exchange is the same as Fig. 3 (a). In4 
(b) a Dy gas with (blue) and without (black) the ther- 
mal exchange is shown. The corresponding condensate 
density is shown with (red) and without (green) thermal 
exchange. 4 (c) Shows a comparison of Dy and Cr with 
a large number of atoms (3700), when the chemical po- 
tentials are nearly the same at zero temperature. The 
densities are normalized to the maximum density at zero 
temperature. 

In 4 (a) the impact of the thermal exchange is not too 
pronounced. It slightly increases the condensate central 
density, which is most obvious at high temperature. This 
is significantly different from the Dy example shown in 4 
(b). There is a noticeable difference between the density 
profiles with and without the thermal exchange. The ef- 
fect is most clearly seen in the condensate by looking at 
the increased central density and reduced extent of the 
red curve (with) versus green (without). At the highest 
temperature, where the validity of HFBP is questionable, 
we see a significant reduction in the condensate density 
and occupation. In fact the condensate repels the ther- 
mal cloud so strongly, that in the center of the trap there 
is a local minimum in the total density. Furthermore, 
the thermal exchange has significantly lowered the con- 
denstate occupation. Thus the inclusion of the thermal 
exchange clearly lowers the critical temperature of the 



dipolar gas. 

In 4 (c) we show that by varying g d and N while hold- 
ing Ng d constant, the HFBP does not give identical re- 
sults as the T = CPE would. For g d = 0.034, N = 300 
(Dy, black) and g d = 0.0028, N = 3700 (Cr, blue) at low 
temperatures the density profiles are slightly different. 
More importantly, at high temperature the profiles are 
very different, this has to do with the strong depletion of 
the condensate. 

In Fig. 5 we show the excitation spectrum of the dipo- 
lar gas (blue circles), a dipolar gas without the thermal 
exchange (red triangle), and the contact gas (black +) 
as a function of temperature. We have also character- 
ized the excitations by their azimuthal symmetry. In a 
2D contact gas, there is a hidden symmetry which makes 
the breathing mode (m=0) have an energy of 2Hco [38]. 
This mode has very little temperature dependence. This 
hidden symmetry is removed by the dipolar interaction, 
however if l z goes to zero, the breathing mode goes to 
2Huj. In the example we have picked (g d — 0.0028 and 
A^=2000), the \m\ = 3 mode is near 2fkv, but this is not 
the breathing mode (to = 0). Rather, the mode just 
below this, is the breathing mode. The contact gas ex- 
citation spectrum (black +) agrees well with Ref. [39]. 
They used g ~ 0.002 which accounts for the differences. 
It is a general feature that the dipolar gas has lower exci- 
tation energies than the contact gas. In Ref. [27], the 3D 
excitation spectrum as a function of temperature showed 
that for an oblate geometry the dipolar excitations are 
lower than the contact system. It is worth pointing, out 
that the dipolar calculations with and without the ther- 
mal exchange, agree well at low temperature. As the 
temperature is increased there is some disagreement be- 
tween the two (blue circle and red triangle), but mostly 
for more highly excitated modes. The calculations with- 
out the exchange become higher in energy. 

The modes that have strong temperature dependence 
are those with higher azimuthal symmetry In contrast, 
the full bodied modes (large central amplitude) have a 
more constant excitation frequency as temperature is var- 
ied. The reason for this is discussed below. 

In Fig. 6 we show the excitation spectrum as a function 
of N for the dipolar (blue circles) at T/T c0 = 0.05. The 
plots on the right are contour plots of the quasi-particle 
modes (u a ) for N — 5000 dipolar atoms. The contour 
lines are in 0.25 increments of the maximum value of u a 
(i.e., 0.75, 0.5, ...,-0.75). Additionally, the negative re- 
gions are shaded blue. The condensate density contours 
are shown as a dashed black lines. The azimuthal sym- 
metries are: (a) m~2 quadruple mode, (b) to=3, (c) m=0, 
breathing mode, (d) m=4, (c) m=l, and (g) m=5. We 
have shown one of the two degenerate modes (±to) at 
each energy when to ^ 0. If the contact system were 
plotted, we would see that excitations are always higher 
in energy for all N shown. 

The modes that are not strongly dependent on num- 
ber or temperature are the modes that have a full bodied 
motion or large amplitude in the middle of the conden- 
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FIG. 3: The total and condensate density of (a) a Cr dipolar gas and (b) contact gas, and the comparison of the total densities 
(c). (a) The total density is in blue (black) for the dipolar (contact) gas, and the condensate density is in red. The temperatures 
from top to bottom are T/T c o — 0.25 (dashed), 0.55 (dotted), 0.75 (dash-dot), and 0.90 (solid), (c) The dipolar gas is more 
dense in the center of the trap at all temperatures. 




FIG. 4: Comparison of density with and without thermal exchange. Both the total and condensate density are shown at four 
temperatures, top to bottom: T/T c o = 0.25 (dashed), 0.55 (dotted), 0.75 (dash-dot), and 0.90 (solid), (a) A Cr gas with (blue) 
and without (black) the thermal exchange; the corresponding condensate density is with (red) and without (green) thermal 
exchange is shown, (b) For a Dy gas, the total density with (blue) and without (black) the thermal exchange; the corresponding 
condensate density with (red) and without (green) thermal exchange is shown. Both the condensate and total density without 
thermal exchange have a lower central density, (c) Shows a comparison of Dy (black) and Cr (blue) with a large number of 
atoms (3700) when the chemical potential are the same at zero temperature. 



sate, see 6 (c) and (e). These modes are mildly impacted 
by the total number of atoms in the system. The modes 
with higher azimuthal symmetries are more like surface 
modes, and are strongly effected by the number or tem- 
perature, see 6 (a), (b), (d), and (f). 

In fact the strong number dependence leads to excita- 
tions crossing paths as N is increased. The m=3, (b), 
mode becomes lower than the breathing mode, (c); and 
m=4, (d), mode becomes lower than the m=l, (e). The 
fact that surface modes move relative to full body excita- 
tions as a function of N is related to the size of the con- 
densate. Both as the atom number is increased, and as 
the temperature is decreased, the size of the condensate 



becomes larger. Naively, one might consider the surface 
of the condensate a ring with a restorating force. We 
consider the surface excitations as displacements of the 
ring from equilibrium, and find an excitation spectrum 
which behaves as m/R where R is radius of the ring [40]. 
Then we look at m over the radius at which the density 
equals 0.25?io(0) (the furthest out contour of the conden- 
sate in Fig. 6), we find a similar number dependence to 
the excitation spectra in Fig. 6. This is too simplistic, 
but gives a physical reason for such behavior in the exci- 
tation spectrum. The important point of this analysis is 
that the size of the condensate grows quickly at small N 
and slows down at larger N. Thus we see a rapid change 
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FIG. 5: The excitation spectrum of the dipolar (blue cir- 
cles), dipolar without exchange (red triangle), and contact 
gas (black +) as a function of temperature for N = 2000. 
The azimuthal symmetry of the excitations are shown next 
to each curve. 

in the excitation spectrum at low N and less so at large 
N. 



VII. CONCLUSION 

This paper studied a trapped q2D dipolar gases at fi- 
nite temperature. We presented the numerical method 
used to solve Hartree Fock Bogoliubov Popov equations 
for a gas with non-local interactions including thermal 
exchange effects. In Fig. 1 we showed condensate frac- 
tion as a function of temperatue for both the Cr and Dy 
examples. The critical temperature of the Dy is greatly 
reduced by the interactions. In Fig. 2 we looked at the 
chemical potential and the total interaction energy as a 
function of temperature for the various terms in Eq. 1. 
We compared the contact example with the Cr examle 
with and without thermal exchange. In Fig. 3 we stud- 
ied the impact of temperature on the density profiles and 



saw that the dipolar gas is more dense in the center of 
the trap. 

Next in Fig. 4 we explored the impact of the thermal 
exchange on both the Cr and Dy examples. We found 
that the strongly interaction Dy example, the thermal ex- 
change strongly reduces the condensate fraction at high 
temperature. In Fig. 5 we studied the impact of tem- 
perature on the excitation spectra. We compared the 
contact example with the Cr examle with and without 
thermal exchange. Finally, in Fig. 6 we looked at the 
excitation spectrum of a Cr g funciton of atom 

number and the quasi-particles, u a , for N = 5000. Here 
we illustrated the strong number or size dependence of 
the high azimuthal symmetries. 

This work has set the stage to study phase coherence 
of a dipolar gas as a function of temperature, as has 
been done for contact gases [39] . We seek to understand 
how the non-local interaction will alter the correlations 
at the BKT transition. We already saw an increased 
phase space density in the middle of the cloud from the 
dipolar interaction. How does the dipolar interaction im- 
pact the phase coherence of the gas? Second, we wish to 
study roton physics in q2D [28, 41, 42]. In q2D, a roton 
can emerge as the field is tilted into the plane of motion, 
and leads to anisotropic density profiles [21]. The pre- 
sented numerical method has been developed to handle 
this configuration. 

This method could be applied to other momentum de- 
pendent interactions, for instance those with renormal- 
ized contact interactions, such as g —> (g + g 2 /k 2 ) [25]. 
Additionally, we could improve the scattering model to 
include more dipolar scattering physics [43]. 
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We present finite temperature analysis of a quasi2D dipolar gas. To do this, we use the Hartree 
Fock Bogoliubov method within the Popov approximation. This formalism is a set of non-local 
equations containing the dipole-dipole interaction and the condensate and thermal correlation func- 
tions, which are solved self-consistently. We detail the numerical method used to implement the 
scheme. We present density profiles for a finite temperature dipolar gas in quasi2D, and compare 
these results to a gas with zero-range interactions. Additionally, we analyze the excitation spectrum 
and study the impact of the thermal exchange. 
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I. INTRODUCTION 

Two hot topics in ultracold physics are two dimen- 
sional (2D) gases and dipolar gases. Reduced dimen- 
sionality enhances the quantum mechanical character of 
the system. In a homogeneous 2D system, quantum 
phase fluctuations are so strong that at finite tempera- 
ture phase coherence cannot be established and conden- 
sation does not occur. There is, however, an intrigu- 
ing phase transition called the Berezinskii-Kostcrlitz- 
Thouless (BKT) transition [1] which occurs when the 
temperature is lowered and there is no longer enough 
thermal energy to unbind vortex and anti-vortex pairs. 
This binding of vortices reduces phase fluctuations so 
that quasi-long range order can form [2]. Interestingly, 
a trapped 2D gas can Bose Einstein condense (BEC) at 
finite temperature [2, 3]. The BKT and BEC transi- 
tion have been studied in trapped ultracold gases, first 
by observations of phase defects [4] , direct observation of 
vortices [5], and then changes in the density profile due to 
the onset of superfluidity [6] . Additionally, universality 
has been observed near the BKT transition [7] . Success- 
ful methods which study such systems are Monte Carlo 
[8, 9], mean field [11, 12] and classical-field methods [13- 
15]. 

In dipolar systems the interactions are non-local and 
the possibility of creating strongly correlated gases are 
tantalizing, especially in quasi-2D (q2D) where zero point 
motion in the axial direction is included. For example, 
some studies have used Monte Carlo methods to study 
phase transitions such as crystallization [16]. This study 
included temperature as a variable and observed melting 
of the dipolar crystal. The BKT transition of a homo- 
geneous dipolar gas has been studied with Monte Carlo 
methods to examine superfluid fraction and excitation 
spectra [17]. Other studies have looked at the structure 
of the dipolar gas and the impact of temperature on the 
phase of the gas [18, 19]. Such studies have focused on 
strongly interacting dipolar systems. However, the bulk 
of the work on q2D gases has been at zero temperature. 
For example, phonon instability [20] and anisotropic su- 
perfluidity [21] have been predicted. 

Recently, chromium (Cr) and dysprosium (Dy) have 



been Bose condensed and exhibited strong dipolar effects 
[22, 23]. Additionally, progress towards a q2D dipolar 
gas has been made with a layered dipolar system for Cr 
atoms in a one dimensional optical lattice [24]. With 
such progress on dipolar gases, direct study of q2D dipo- 
lar systems has begun. An important question is how 
will dipolar interactions impact the quantum behavior 
of a q2D gas? The long range nature of the dipolar in- 
teraction may lead to interesting physics for the thermal 
system, especially relating the BKT transition and phase 
coherence. Additionally, little work studying the impact 
of temperature on q2D dipolar gases has been conducted 
for reasonable experimental parameters. 

With an eye toward this unexplored physics and aiding 
experiments, we study the finite temperature physics of 
q2D trapped dipolar gases. We use the standard beyond 
mean-field method: the Hartree Fock Bogoliubov within 
the Popov approximation (HFBP) [25]. The HFBP has 
been successful in studying ultracold atoms, for example 
it has been used to explain the temperature dependence 
of collective excitations in a 3D BEC [26]. In the case 
of 3D dipolar gases, Rcf. [27] used the HFBP to study 
temperature effects on the biconcave structure dipolar 
gases can have [28]. This study neglected the thermal 
exchange. Other recent work used mean field methods 
to study the stability of finite temperature dipolar gases 
[29]. 

The aim of this paper is two fold. First, we develop 
a numerical method to solve the non-local HFBP. This 
method includes the non-local interaction and exchange 
effects. We show in detail how the interaction and ther- 
mal exchange effects arc included; it relies on a parallel 
implementation of the method. The method can also 
be applied to dipolar fermions. Second, we compare the 
contact gas and the dipolar gas with the HFBP. In the 
comparison of the two gases, we look at the condensate 
number as a function of temperature. We also articulate 
the role of thermal exchange in the gas by solving the sys- 
tem with and without thermal exchange. We compare 
the density profiles of the gas at various temperatures. 
We look at the excitation spectrum as a function of both 
temperature and total particle number. We also classify 
the lowest excitation modes. To the author's knowledge, 
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these results are the first for a q2D trapped-finite tem- 
perature dipolar gas with thermal exchange effects. 

II. EQUATIONS OF MOTION 

We study a quasi-2D dipolar system at finite tempera- 
ture. To do this, we will employ the Hartree Fock Bogoli- 
ubov method within the Popov approximation with non- 
local interactions [25, 30]. The HFB breaks the second 
quantized bosonic field operator into a condensate and 
non-condensate (thermal) component: ^ = [^/Wo(f) {x) + 



9{x)] where we have replaced So — > \/N a . Here x repre- 
sents all required coordinates, for the 2D case it will be 
p. We will use the Bogoliubov transformation for the 
thermal part: 6(x) ~ J2 a i u »( x )^ae^' lu ' a ' t ~ v^(x)d* a e lulat ] 
where d (d*) is the bosonic annihilation (creation) oper- 
ator for a quasi-particle (hole). They obey the bosonic 
commutation relations: [a a ,aV\ = <5 Q( g and [a,,,^] = 
[d*,dg] = 0. To derive the equations of motion we start 
with the Heisenberg equation of motion with a non-local 
interaction [30]. We find for our system the modified 
non-local Gross Pitaevskii equation (GPE) is 



H<p (x) = (h + J dx'n(x')V^j <f) (x) + J dx'lhix,^) - m(x,x')]V(j> (x'), (1) 

where Hq is the kinetic energy and trapping potential. The interaction potential, V(x — x'), has been written as V 
for simplicity. The total density, n(x) = uq(x) + fi{x) is made up of condensate density and thermal density. The 
non-local correlation functions are: no(x,x') = N <f)Q(x)(f)o{x') and n (x) = n (x 7 x); n(x,x') = (9* (x)9(x')) is the 
thermal correlation function and the local thermal density is h{x) = hix^x). rh{x,x') = (9(x)9(x')} is the anomalous 
thermal correlation function, and in the Popov approximation is neglected. The quasi-particle wavefunctions are 
required to construct the thermal and anomalous correlation functions. At equilibrium, these wavefunctions can be 
found from a non-local Bogoliubov-de Gennes equation: 



+ J dx'm(x, x')Vv a {x') — J dx'm (x,x')Vv a (x') 



(2) 



-u n v n 



+ 



J dx'rh*(x, x')Vu a (x') — J dx'm^(x,x')Vu a (x'). 



The total number of atoms is N = N + N and N — 
J dxn(x). The normalization of the quasi-particle wave- 
function is 1 = J dx(\u a \ 2 — \v a \ 2 ). The anomalous con- 



densate correlation function is m (x,x') = N <j)(x)(j)(x'). 
The thermal correlation functions are defined in terms of 
the quasi-particles: 



h(x,x') = J2K( X > a ( X ) + v a (x')v* a (x)}N^ + v a (x')v* a (x) (3) 

a 

m(x,x') =J2K(x')v* a (x)+v* a (x')u a (x)]N? e + u a (x')v* a (x), (4) 

i 

where Ng e = (Ze hUa / kT — 1) _1 = (a* d a ) where k is cation is made. 
Boltzmann's constant, T is the temperature, and Z = 
1 + 1 /No . This relation is simply the Bosonic occupation 
of thermal modes. These equations are consistent with 
Refs. [27], [30], and [31] when the appropriate simplifi- 
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To make the discussion of Eqs. (1) and (2) simpler, we 
introduce a compact notation: 



is: 



{H + D + X + M)fa= vufo, 
(H Q - p, + D + X)u a + Mv a 
(H - p. + D + X)v a + Mu a 



(5) 
(6) 



-i0 a V o 



The first term is the standard kinetic and potential 
energy, Hq. For the remaining interaction terms, we 
use the capital letters to indicate the inclusion of the 
x' integral, so the direct interaction is D = Dq + D 
(condensate and thermal), and for example D (x) = 
J dx'[no(x')]V(x — x'). The exchange interaction is 
X = Xq + X (condensate and thermal), and for example 
X(x,x') — J dx'[h(x,x')]V(x — x'). Finally, the anoma- 
lous correlation function is M = M—M , and for example 
M(x,x')= J dx'[rh(x,x') - m (x,x')]V(x - x'). We will 
use the Popov approximation, which assumes M = 0, 
leading to M — — A if (f>o is real. This approximation 
keeps the spectrum gapless. 

We project out the condensate mode from the quasi- 
particles, as was done in Refs. [32, 33] which uses 
Q = 1 — |0)(0| and (a;|0) = (j>o(x). The projection op- 
erator is applied to all operators that are not in H GP = 
Hq — /i + D + X, such as Mq and Xq, and for exam- 
ple, Mq — > QM Q. This projection method keeps the 
spectrum gapless. 



III. THE INTERACTION 

We are interested in the q2D dipolar system assuming 
no axis of symmetry, so we take the dipole moment to 
be d = d(zcos(a) + xsm(a)), and we assume harmonic 
confinement and that ui z 3> lo p where the uii are the 
trapping frequencies ( u> x = ui y = u> p ). This allows us 
to assume only one transverse mode is occupied in the 
z direction. We can then evaluate the interaction by 
factoring the wavefunction as ij)(r) — xo(z)ip(p), where 
tp is either (f> , u a , or u a and xo is eigenstate of H in 
the z direction. We are then able to integrate out z, and 
obtain an effective interaction. We will also consider a 
contact interaction. The interaction we use is: 



V(p)=gS(p)+g d V dd (p), 



(7) 



where g = \Z8nh 2 a s /ml z is the strength of the contact 
interaction, a s is s-wave the scattering length (a s -C l z ), 
l z = yjh/muj z is the axial harmonic oscillator length, 
gd = d 2 j\phxl z and d is the dipole moment of the par- 
ticles. This interaction leads to the dipolar length scale: 
Id = md 2 /h 2 . In practice, we rescale the equations 
of motion into oscillator units, where t he unit s of en- 
ergy and length are: Hco p and l p = yjh/muj p . After 
this rescaling, the dipolar interaction strength becomes 
ga = ld/{lzV%n)- Finally, the dipolar interaction in q2D 



V dd (p-p)= J dzxl(z) J dz' X 2 (z')V! D (r-n%) 

V2 d D (r)=d 2 (l-3(d.ff)/r*. 

The interaction is evaluated in momentum space. As an 
example, the direct interaction is: 



D(p) = J- 1 \n(k)V (fc) 

T~r/T\ 47T „ / kl z \ 

Vik) = 9d Y F [f2)> 



(9) 



where T is the Fourier transform operator and n(k) = 
F[n(p)\. V k is the momentum representation of the in- 
teraction. The function F(k) is the fc-space dipolar in- 
teraction for the q2D geometry It has two contribu- 
tions coming from polarization perpendicular to and in 
the direction of the dipole tilt, F(k) — cos 2 (a)F±(k) + 
sin 2 (a)F||(fc) where a is the angle between z and the 
polarization vector d. Its contributions arc F\\(k) = 

— 1 + 3v / 7r(fc 2 /fc)e fe erfc(fc), where k d is the wave vec- 
tor along the polarization direction in the x-y plane, 
erfc is the complementary error function and F±(k) = 
2 - 3V^ce fe2 erfc(fc) [20, 21, 34]. 

To evaluate the exchange interaction, we use a simple, 
yet memory intensive, solution: we explicitly construct 
the interaction on the same space-space grid as the cor- 
relation functions. We put n(p,(7) on a grid specified by 
the indices s and t, representing p and p . Each grid has 
n s spatial points, and a correlation function has n 2 grid 
points. 

To evaluate the interaction, we work in momentum 
space. The interaction on the space-space grid is: 

V st = Wj k V k W kt . (10) 

where W ks is the operator that transforms between space 
(s) and k-space (k) representation via the spectral ba- 
sis set [35] . The operator takes the place of the Fourier 
operator in our algorithm. An important aspect of the 
method is that it regularizes the interaction and avoids 
the logarithmic divergence as p — > ff encountered by di- 
rectly evaluating Eq. (8) because W kt involves a project- 
ing onto the Gauss-Hermite basis set [35]. Constructing 
the interaction involves the most costly step, requiring 
rig operations per node. This is performed only once for 
a set of interaction parameters (l z and a) and basis set 
size. 

We have constructed a parallel implementation of the 
method. We distribute the s index (p grid) across the 
computing nodes (p); we will denote this distributed in- 
dex as s p . This way, each core handles only a fraction 
of the correlation functions and interaction. In the fu- 
ture, we will construct the interaction in the partial wave 
expansion: V(p) = ^2 m V m (p 7 (f>). This expansion will 
be necessary to handle an interaction which depends on 
quantum numbers such as n z (confinement in z) or spin. 
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IV. NUMERICAL IMPLEMENTATION 

We solve Eqs. (5) and (6) by expanding the wavefunc- 
tions on the basis that diagonalizes Hq. We will focus 
on the harmonically trapped case where H^Xa = ^aXai 
I dpx* a (p)Xp(p) = <W, and e a = fuv p (m x + m y + 1) 
where m a is an integer. In a basis set, the evaluation 
of the direct interaction term, D a p, is straight forward: 
I dpXa(p)[J dx'V(p- p , )n(p , )]xp(p), where the quantity 
in the square brackets is the effective potential from the 
interaction. The more complicated terms are the non- 
local exchange terms, such as the thermal exchange: X a p 
which is / dpx* a (p) J df?n(p, ff )V(p- ?) X p{P). 

The non-local exchange term in the spectral basis is: 
X aa = U^ s V st n st Uta, where U SIJ is the transformation 
between the spectral basis (a) and space basis (s). The 
numerical procedure to evaluate this on node p is: (1) 
Y Sp t = V Sp tn Spt , multiply the interaction and correlation 
function; (2) h Sp<7 — Y Spt U tcr , project the t spatial ba- 
sis onto the a spectral basis; (3) X^ a = J2 Sp UT p a^s p a, 
project s p spatial basis onto the a spectral basis; and (4) 
Xao- = Sp^ao-i collect and sum. In comparison, the 
direct terms are straight forward and evaluated in mo- 
mentum space at each step. $ s = W s kVkWktnt, where 
& s is the effective potential. Then we must project onto 
the basis: D aa = £/J s 4> s £7 s(7 , and this is done in parallel. 

The procedure to find the full solution is: set temper- 
ature, T, and iVo and pick a targeted value for the total 
number of particles in the system, N tar get, and set h = 
and m — 0: 

1. Solve Eq. (1) for <p and p (in basis set). 

2. Construct condensate exchange term, X , with 
n n (p, /f) and put in basis set. 

3. Solve Eq. (2) for u a and v a . 

4. Construct n(p, ff) and m(p,p'). 

5. Use the semi-classical LDA to supplement the ther- 
mal tail (see below), which gives N and therefore 
N = N + N. 

6. Construct thermal and anomalous exchange terms, 
X and M (if not using Popov approximation), and 
put in basis set representation. 

7. Adjust A to get the desired N tar get of atoms. We 
require N to be within 1% of A 'target ■ 

8. Go back to step 1 until self-consistency is reached. 
We converge the number of thermal atoms, so that 
\hi — fij_i| < 5 • 10~ 5 , when we are near enough 

^target • 

With the numerical method in hand, we are ready to 
proceed to the physical examples and the results. 

The semi-classical, local density approximation (LDA) 
is very important for large p and high momentum states 
on the grid. At high temperature with a manageable 



grid size (much bigger than the condensate), we find that 
there can be a appreciable number of atoms outside of 
quadrature grid. For the large p contributions, the den- 
sity is low enough that the analytic, non-interacting so- 
lution can be used to account for these off-grid particles. 
We integrate the analytic solution over a large region out- 
side of the grid. We have tested that this does not impact 
the convergence of the end result. We include the impact 
of higher transverse modes in the same analytic manner 
on and off the grid. We also assume that these particles 
in the higher z modes are non-interacting. This is rea- 
sonable, especially because symmetry prevents the first 
excited transverse mode from colliding with the conden- 
sate with ui — symmetry. Additionally, a noticable 
number of thermal particles are in momentum states be- 
yond calculated quantum spectrum [36, 37]. To include 
these thermal atoms we use the method from Ref. [27] . 



V. PHYSICAL EXAMPLES 

We will study Cr and Dy atoms. These atoms have 
been Bose condensed and shown to exhibit strong dipolar 
effects [22, 23]. The first example we choose is: 2 • 10 3 Cr 
atoms with Id = 45a where ao is the Bohr radius. This 
gives gd = 0.0028 and g = in trap units. Addition- 
ally, we pick l z /l p = 0.1, this is from using (oj p ,uj z )=2ir 
(16,1600) Hz. We have found that our numerical proce- 
dure using 465 basis states with an energy cut of 30hto p 
converges. To compare this to a contact gas, we match 
the chemical potential at zero temperature by varying g. 
We will only consider a = here, that is perpendicular 
polarization where the dipolar interaction is isotropic and 
repulsive. We also study Dy with a trapping potential of 
(<jJ p ,u z )=2-k (10,1000) Hz. With l d = 400a , this gives 
gd = 0.034 in trap units. This is strongly interacting 
and to maintain the Kohn mode (center of mass slosh 
mode), we have used a low atom number, and for our 
example we use 300 atoms [22]. It is important to note 
that the Cr example is very similar to recent experiments 
[24]. In that work, they were able to vary l z /l p between 
0.1 and 0.17 and have up to 2000 atoms in a layer [24]. 
However, this is a layered system and interlayer dipolar 
interactions are important. 

In contrast to a homogenous system, a trapped ideal 
2D gas can have Bose condensation [3]. In this case, the 
critical temperature is at T c0 /Huj = V&N /n <~ 0.78v^V. 
We will use this temperature to rescale our findings so 
that, to first order, we remove the number dependence 
of the results. This is important because, as we vary 
the temperature, we do not have exactly the same num- 
ber of particles between every calculation. In the ther- 
modynamic limit, the population of the condensate is 
N /N - 1 - (T/T c0 ) 2 [3]. 

For reference, we give several specific numerical exam- 
ples in Table I. For two temperatures, we report p, N /N, 
N, and the Kohn mode energy (should be lfiw p ). For 
the Dy (g d = 0.034) example, at T/T c0 = 0.5, the Kohn 
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TABLE I: Numerical examples of the HFBP for the q2D dipo- 
lar gas and contact gas. 





TIT n 




Nn IN 
iVQ/iV 


N 


rC c\ n n A/f r\c\ p I ni i) 

J.YW11I1 1V1UU.L/ 1 LUJ 


g d =0.0028 


0.05 
0.5 


3.822 
3.457 


0.994 
0.658 


2004 
1994 


0.9992 
0.995 


g d =0.0028 
X — 


0.05 
0.5 


3.814 
3.341 


0.994 
0.660 


2004 
1994 


0.9996 
0.997 


3=0.021 


0.05 
0.5 


3.829 
3.541 


0.993 
0.662 


2004 
1993 


0.9991 
0.9995 


3d =0.034 


0.05 
0.5 


5.124 
4.919 


0.983 
0.525 


300 
298 


0.991 
0.976 


5d=0.034 
X = 


0.05 
0.5 


5.094 
4.524 


0.983 
0.5500 


300 
303 


0.996 
0.987 








FIG. 1: Nq/N as a function of Temperature for the Cr dipolar 
gas (thick black line with circles), the contact gas (green dash 
dot), and for 2000 ideal trapped particles (solid black line). 
Additionally, a gas of 300 Dy particles (dashed red with +) 
and 300 ideal particles are shown (red open squares). 

mode has noticably deviated from 1 by about 2.5%, this 
is about the worst the convergence of this mode gets as a 
function of temperature. For both dipolar examples, we 
show the results for when thermal exchange is negelected 
(X = 0). The lower temperature solutions are similar to 
the GPE solutions at T = 0. 



VI. PROPERTIES OF A Q2D DIPOLAR GAS AT 
FINITE TEMPERATURE 

In this section we study the properties of a q2D dipolar 
gas at finite temperature, and compare it to a contact gas 
and a dipolar system without the thermal exchange. We 
look at several different aspects of these systems. First, 
we look at the condensate number as a function of tem- 
perature. Then, we compare the interaction contribu- 
tions to the total energy. Third, we compare the density 




FIG. 2: The chemical potential for three systems: dipolar 
(black), contact(dashed red), and dipolar without exchange 
(dash dot blue). The interaction terms for each system in 
the Hamiltonian as a function of temperature. The contribut- 
ing terms are the direct condensate ({0o|-Do|</>o)/a*, triangles), 
the direct thermal ({(f)o\D\(f>o) / n, circles), and the thermal ex- 
change {{4> \X /\4>q) / ^l, squares) interaction. 



profiles of the contact and dipolar gas, at various tem- 
peratures. Fourth, we look at the impact of including the 
thermal exchange on the density profile of the gas. Fi- 
nally, we look at the excitation spectrum as a function of 
temperature and total particle number. We also classify 
the lowest excitation modes. 

Figure 1 shows the condensate fraction for the Cr dipo- 
lar gas (black line with circles), the contact gas (green 
dash dot line), for 2000 ideal particles (solid black line), 
300 Dy particles (dashed red), and 300 ideal particles (red 
open squares). The condensate fraction for the dipolar 
gas and contact gas are shifted down from the ideal gas 
(black line), with the dipolar gas being slightly lower than 
the contact gas. For the Dy example, it is clear that the 
interaction strongly depletes the condensate mode. 

In Fig. 2 we look at the chemical potential and its 
contributions from Eq. 1, as a function of temperature. 
The examples shown are dipolar (black), contact (red), 
and dipolar without thermal exchange (blue). We show 
the direct condensate interaction ((0o|-Do|0o)/M' trian- 
gle), the direct thermal interaction (((j>o\D\(j>o) / /i, circle), 
and the thermal exchange interaction ((c/)o\X /\<t>o)/ /i, 
square). For the contact interaction, the thermal ex- 
change is equal to thermal direct interaction, so we only 
show one. For each example, the remaining contribution 
to the chemical potential is from Hq or the potential and 
kinetic energy contributions. 

For the dipolar gas, the thermal exchange term is 
about half as strong the direct contribution. The im- 
portance of the interaction with the thermal particles is 
more important as the temperature is increased. Com- 
paring the dipolar calculation with and without the ther- 
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mal exchange, we see that the chemical potentials agree 
reasonably well, the full calculation has a larger /i at 
high temperature. Looking at the interaction contribu- 
tions, we see that the direct condensate (black triangle) 
and thermal interaction (black circles) make up a smaller 
portion of the chemical potential for the full calculation 
than the one with X = (blue triangle and open circles) . 

In Fig. 3 we show both the total and condensate den- 
sity for (a) a Cr system of 2000 particles and an analogous 
(b) contact system (near equal chemical potentials). We 
then compare the density profiles in (c). The total den- 
sity is in blue (black) for the dipolar (contact) systems, 
and the condensate density is shown in red. The temper- 
atures from top to bottom are T/T c0 — 0.25 (dashed), 
0.55 (dotted), 0.75 (dash-dot), and 0.90 (solid). For both 
the contact and dipolar gas at large p, the thermal den- 
sity behaves as the analytic solution predicts. 

An important point is that the dipolar gas has a higher 
density in the middle of the trap than the contact gas. It 
is hard to see in the figure, but the contact condensate 
atoms have been shifted to the shoulders of the conden- 
sate or near the trap edge. This has implications for the 
temperature at which the superfiuid phase transition will 
occur in a dipolar gas. 

In Fig. 4 we compare both the total and condensate 
density with and without the thermal exchange. We 
look at four temperatures, from top to bottom they are: 
T/T c0 = 0.25 (dashed), 0.55 (dotted), 0.75 (dash-dot), 
and 0.90 (solid). In 4 (a), a Cr gas with (blue) and with- 
out (black) the thermal exchange is shown. Additionally, 
the corresponding condensate density is with (red) and 
without (green) thermal exchange. This example with 
the thermal exchange is the same as Fig. 3 (a). In4 
(b) a Dy gas with (blue) and without (black) the ther- 
mal exchange is shown. The corresponding condensate 
density is shown with (red) and without (green) thermal 
exchange. 4 (c) Shows a comparison of Dy and Cr with 
a large number of atoms (3700), when the chemical po- 
tentials are nearly the same at zero temperature. The 
densities are normalized to the maximum density at zero 
temperature. 

In 4 (a) the impact of the thermal exchange is not too 
pronounced. It slightly increases the condensate central 
density, which is most obvious at high temperature. This 
is significantly different from the Dy example shown in 4 
(b). There is a noticeable difference between the density 
profiles with and without the thermal exchange. The ef- 
fect is most clearly seen in the condensate by looking at 
the increased central density and reduced extent of the 
red curve (with) versus green (without). At the highest 
temperature, where the validity of HFBP is questionable, 
we see a significant reduction in the condensate density 
and occupation. In fact the condensate repels the ther- 
mal cloud so strongly, that in the center of the trap there 
is a local minimum in the total density. Furthermore, 
the thermal exchange has significantly lowered the con- 
denstate occupation. Thus the inclusion of the thermal 
exchange clearly lowers the critical temperature of the 



dipolar gas. 

In 4 (c) we show that by varying g d and N while hold- 
ing Ng d constant, the HFBP does not give identical re- 
sults as the T = CPE would. For g d = 0.034, N = 300 
(Dy, black) and g d = 0.0028, N = 3700 (Cr, blue) at low 
temperatures the density profiles are slightly different. 
More importantly, at high temperature the profiles are 
very different, this has to do with the strong depletion of 
the condensate. 

In Fig. 5 we show the excitation spectrum of the dipo- 
lar gas (blue circles), a dipolar gas without the thermal 
exchange (red triangle), and the contact gas (black +) 
as a function of temperature. We have also character- 
ized the excitations by their azimuthal symmetry. In a 
2D contact gas, there is a hidden symmetry which makes 
the breathing mode (m=0) have an energy of 2Hco [38]. 
This mode has very little temperature dependence. This 
hidden symmetry is removed by the dipolar interaction, 
however if l z goes to zero, the breathing mode goes to 
2Huj. In the example we have picked (g d — 0.0028 and 
A^=2000), the \m\ = 3 mode is near 2fkv, but this is not 
the breathing mode (to = 0). Rather, the mode just 
below this, is the breathing mode. The contact gas ex- 
citation spectrum (black +) agrees well with Ref. [39]. 
They used g ~ 0.002 which accounts for the differences. 
It is a general feature that the dipolar gas has lower exci- 
tation energies than the contact gas. In Ref. [27], the 3D 
excitation spectrum as a function of temperature showed 
that for an oblate geometry the dipolar excitations are 
lower than the contact system. It is worth pointing, out 
that the dipolar calculations with and without the ther- 
mal exchange, agree well at low temperature. As the 
temperature is increased there is some disagreement be- 
tween the two (blue circle and red triangle), but mostly 
for more highly excitated modes. The calculations with- 
out the exchange become higher in energy. 

The modes that have strong temperature dependence 
are those with higher azimuthal symmetry In contrast, 
the full bodied modes (large central amplitude) have a 
more constant excitation frequency as temperature is var- 
ied. The reason for this is discussed below. 

In Fig. 6 we show the excitation spectrum as a function 
of N for the dipolar (blue circles) at T/T c0 = 0.05. The 
plots on the right are contour plots of the quasi-particle 
modes (u a ) for N — 5000 dipolar atoms. The contour 
lines are in 0.25 increments of the maximum value of u a 
(i.e., 0.75, 0.5, ...,-0.75). Additionally, the negative re- 
gions are shaded blue. The condensate density contours 
are shown as a dashed black lines. The azimuthal sym- 
metries are: (a) m~2 quadruple mode, (b) to=3, (c) m=0, 
breathing mode, (d) m=4, (c) m=l, and (g) m=5. We 
have shown one of the two degenerate modes (±to) at 
each energy when to ^ 0. If the contact system were 
plotted, we would see that excitations are always higher 
in energy for all N shown. 

The modes that are not strongly dependent on num- 
ber or temperature are the modes that have a full bodied 
motion or large amplitude in the middle of the conden- 
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FIG. 3: The total and condensate density of (a) a Cr dipolar gas and (b) contact gas, and the comparison of the total densities 
(c). (a) The total density is in blue (black) for the dipolar (contact) gas, and the condensate density is in red. The temperatures 
from top to bottom are T/T c o — 0.25 (dashed), 0.55 (dotted), 0.75 (dash-dot), and 0.90 (solid), (c) The dipolar gas is more 
dense in the center of the trap at all temperatures. 




FIG. 4: Comparison of density with and without thermal exchange. Both the total and condensate density are shown at four 
temperatures, top to bottom: T/T c o = 0.25 (dashed), 0.55 (dotted), 0.75 (dash-dot), and 0.90 (solid), (a) A Cr gas with (blue) 
and without (black) the thermal exchange; the corresponding condensate density is with (red) and without (green) thermal 
exchange is shown, (b) For a Dy gas, the total density with (blue) and without (black) the thermal exchange; the corresponding 
condensate density with (red) and without (green) thermal exchange is shown. Both the condensate and total density without 
thermal exchange have a lower central density, (c) Shows a comparison of Dy (black) and Cr (blue) with a large number of 
atoms (3700) when the chemical potential are the same at zero temperature. 



sate, see 6 (c) and (e). These modes are mildly impacted 
by the total number of atoms in the system. The modes 
with higher azimuthal symmetries are more like surface 
modes, and are strongly effected by the number or tem- 
perature, see 6 (a), (b), (d), and (f). 

In fact the strong number dependence leads to excita- 
tions crossing paths as N is increased. The m=3, (b), 
mode becomes lower than the breathing mode, (c); and 
m=4, (d), mode becomes lower than the m=l, (e). The 
fact that surface modes move relative to full body excita- 
tions as a function of N is related to the size of the con- 
densate. Both as the atom number is increased, and as 
the temperature is decreased, the size of the condensate 



becomes larger. Naively, one might consider the surface 
of the condensate a ring with a restorating force. We 
consider the surface excitations as displacements of the 
ring from equilibrium, and find an excitation spectrum 
which behaves as m/R where R is radius of the ring [40]. 
Then we look at m over the radius at which the density 
equals 0.25?io(0) (the furthest out contour of the conden- 
sate in Fig. 6), we find a similar number dependence to 
the excitation spectra in Fig. 6. This is too simplistic, 
but gives a physical reason for such behavior in the exci- 
tation spectrum. The important point of this analysis is 
that the size of the condensate grows quickly at small N 
and slows down at larger N. Thus we see a rapid change 
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FIG. 5: The excitation spectrum of the dipolar (blue cir- 
cles), dipolar without exchange (red triangle), and contact 
gas (black +) as a function of temperature for N = 2000. 
The azimuthal symmetry of the excitations are shown next 
to each curve. 

in the excitation spectrum at low N and less so at large 
N. 



VII. CONCLUSION 

This paper studied a trapped q2D dipolar gases at fi- 
nite temperature. We presented the numerical method 
used to solve Hartree Fock Bogoliubov Popov equations 
for a gas with non-local interactions including thermal 
exchange effects. In Fig. 1 we showed condensate frac- 
tion as a function of temperatue for both the Cr and Dy 
examples. The critical temperature of the Dy is greatly 
reduced by the interactions. In Fig. 2 we looked at the 
chemical potential and the total interaction energy as a 
function of temperature for the various terms in Eq. 1. 
We compared the contact example with the Cr examle 
with and without thermal exchange. In Fig. 3 we stud- 
ied the impact of temperature on the density profiles and 



saw that the dipolar gas is more dense in the center of 
the trap. 

Next in Fig. 4 we explored the impact of the thermal 
exchange on both the Cr and Dy examples. We found 
that the strongly interaction Dy example, the thermal ex- 
change strongly reduces the condensate fraction at high 
temperature. In Fig. 5 we studied the impact of tem- 
perature on the excitation spectra. We compared the 
contact example with the Cr examle with and without 
thermal exchange. Finally, in Fig. 6 we looked at the 
excitation spectrum of a Cr g funciton of atom 

number and the quasi-particles, u a , for N = 5000. Here 
we illustrated the strong number or size dependence of 
the high azimuthal symmetries. 

This work has set the stage to study phase coherence 
of a dipolar gas as a function of temperature, as has 
been done for contact gases [39] . We seek to understand 
how the non-local interaction will alter the correlations 
at the BKT transition. We already saw an increased 
phase space density in the middle of the cloud from the 
dipolar interaction. How does the dipolar interaction im- 
pact the phase coherence of the gas? Second, we wish to 
study roton physics in q2D [28, 41, 42]. In q2D, a roton 
can emerge as the field is tilted into the plane of motion, 
and leads to anisotropic density profiles [21]. The pre- 
sented numerical method has been developed to handle 
this configuration. 

This method could be applied to other momentum de- 
pendent interactions, for instance those with renormal- 
ized contact interactions, such as g —> (g + g 2 /k 2 ) [25]. 
Additionally, we could improve the scattering model to 
include more dipolar scattering physics [43]. 
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